Spatial Ecology and Macroecology

Practical - Week 2

Florencia Grattarola, Friederike Woelke & Gabriel Ortega

(Department of Spatial Sciences)

2022-10-10

What are we going to see today?

We will explore the potential drivers of species distribution.

  1. Overview of data sources
  2. Get data
    1. Process raster
    2. Visualise data
  3. Correlations between predictors

What are we going to see today?

We will introduce three new packages.

  1. terra: tools to manage and analyze rasters and vectors

  2. geodata: interface to multiple sources of spatial data

  3. MODIStsp: pipelines for processing time series from MODIS data

worldclim.org

WorldClim Is a database of high spatial resolution global weather and climate data. You can download gridded weather and climate data for historical (near current) and future conditions.

chelsa-climate.org

CHELSA (Climatologies at high resolution for the earth’s land surface areas) is a very high resolution (30 arc sec, ~1km) global downscaled climate data set. It includes climate layers for various time periods and variables, ranging from the Last Glacial Maximum, to the present, to several future scenarios.

nature.com/articles/s41597-020-00599-8

A global, spatially explicit characterization of 47 terrestrial habitat types, as defined in the International Union for Conservation of Nature (IUCN) habitat classification scheme.

Habitatmapping: https://github.com/Martin-Jung/Habitatmapping

ecoregions.appspot.com

The RESOLVE Ecoregions dataset, updated in 2017, offers a depiction of the 846 terrestrial ecoregions that represent our living planet.


### land.copernicus.eu/pan-european/corine-land-cover
The CORINE Land Cover (CLC) consists of an inventory of land cover in 44 classes (from Europe). The inventory was initiated in 1985. Updates have been produced in 2000, 2006, 2012, and 2018.

modis.gsfc.nasa.gov

Terra MODIS and Aqua MODIS satellites are viewing the entire Earth’s surface every 1 to 2 days, acquiring data in 36 spectral bands or groups of wavelengths.

MODIStsp: https://github.com/ropensci/MODIStsp/

sedac.ciesin.columbia.edu

SEDAC, the Socioeconomic Data and Applications Center, focuses on human interactions in the environment and seeks to develop and operate applications that support the integration of socioeconomic and earth science data.

EXCERCISE Drivers

Drivers

In today’s practical, we will focus on Indonesia and assess five different layers:

  • Annual mean temperature (WorldClim BIO1)
  • Precipitation seasonality (WorldClim BIO15)
  • Elevation (SRTM)
  • Land-cover (MODIS Terra MCD12Q1)
  • Human footprint index (SEDAC)

2. Setup your work environment: libraries

2.0 Install pacman

pacman is a handy package to install, load, and unload libraries. If you attempt to load a library that is not installed, pacman will try to install it automatically.

install.packages("pacman")

2.1 Load your libraries

We will need tidyverse, geodata, terra, MODIStsp, and sf

pacman::p_load(
  rnaturalearth, # useful base maps
  geodata, # Get data from multiple sources
  terra, # Work with rasters and vectorial spatial data
  MODIStsp, # Get data from MODIS
  sf,  # Manage spatial data as simple features dataframes
  tmap, # Cool maps
  tictoc # Benchmarking
)

pacman::p_load(tidyverse) # Data wrangling... and.. why not?

2.2 Create project variables

We will need:

  1. The polygon and spatial extent of Indonesia.

  2. A planar projection of the area.

  3. The taxa we are going to analyze.

2.2 Create project variables

Get a polygon of Indonesia.

country <- ne_countries(
  country = "indonesia",
  returnclass = "sf"
) # We want an object of class sf

2.2 Create project variables

Print the extent of the country plus 1 degree. If you work in a small area, you cannot get rid of decimal values in coordinates. However, working in a large area you can round the values to simplify your life!.

round(st_bbox(country),0) # Check the extent of the country (round values)
xmin ymin xmax ymax 
  95  -10  141    5 
round(st_bbox(country),0) + c(1,-1,1,1) #This way we add space around the country
xmin ymin xmax ymax 
  96  -11  142    6

Now visit https://projectionwizard.org/ and use the second set of coordinates above to find the recommended projection for your work today.

2.2 Create project variables

Now, let’s create an object containing the project CRS. We will not transform our spatial objects to this CRS immediately, but it will be useful later.

project_crs <- 'PROJCS["ProjWiz_Custom_Cylindrical_Equal_Area",
 GEOGCS["GCS_WGS_1984",
  DATUM["D_WGS_1984",
   SPHEROID["WGS_1984",6378137.0,298.257223563]],
  PRIMEM["Greenwich",0.0],
  UNIT["Degree",0.0174532925199433]],
 PROJECTION["Cylindrical_Equal_Area"],
 PARAMETER["False_Easting",0.0],
 PARAMETER["False_Northing",0.0],
 PARAMETER["Central_Meridian",119],
 PARAMETER["Standard_Parallel_1",0],
 UNIT["Meter",1.0]]'

2.2 Create project variables

Create an extent around Indonesia.

extent <- round(st_bbox(country),0) + c(1,-1,1,1) # Bounding box plus 1 degree

2.2 Create project variables

Set the path were you are going to save data. This is required to use the package geodata. From now on, every time you try to get a dataset with geodata functions, the package will check if the data already exists in your data folder and load it from there.

geodata_path("data/")

3. Get raster data

Recommendations

When working with rasters, there are three basic procedures we can recommend you:

  1. Crop: Cut the rasters to the geographic extent you need. This is a square cut.

  2. Reproject: Make sure that all your rasters are on the same projection. Be aware that reprojecting raster is not the best. Try to reproject your vector data to match the raster.

  3. Mask: Set to NA all the cells outside your polygon(s) of interest.

3. Get data: climatic

Let’s download the WorldClim data (version 2.1 climate data for 1970-2000):

  • Annual mean temperature (bio1)
  • Annual precipitation (bio15)


You can download these data from https://www.worldclim.org/data/worldclim21.html
or use the package geodata.

2. Get data: climatic

We will use the function worldclim_country() from geodata.

tavg <- worldclim_country(
  var = "tavg",
  country = "IDN" # IDN is the three letters ISO code for Indonesia.
)

tavg <- crop(tavg, extent)

|---------|---------|---------|---------|
=========================================
                                          
tavg
class       : SpatRaster 
dimensions  : 2040, 5460, 12  (nrow, ncol, nlyr)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : 96, 141.5, -11, 6  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : spat_Bj682IhX8aXcA0V_59366.tif 
names       : IDN_w~avg_1, IDN_w~avg_2, IDN_w~avg_3, IDN_w~avg_4, IDN_w~avg_5, IDN_w~avg_6, ... 
min values  :         3.9,         3.7,         3.8,         3.8,         3.8,         2.9, ... 
max values  :        30.2,        29.9,        29.6,        29.3,        29.7,        29.7, ...

2. Get data: climatic

Now get the precipitation data from Worldclim.

prec <- worldclim_country(
  var = "prec",
  country = "IDN"
)

prec <- crop(prec, extent)

|---------|---------|---------|---------|
=========================================
                                          
prec
class       : SpatRaster 
dimensions  : 2040, 5460, 12  (nrow, ncol, nlyr)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : 96, 141.5, -11, 6  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : spat_qeaawbMRAtGAzfP_59366.tif 
names       : IDN_w~rec_1, IDN_w~rec_2, IDN_w~rec_3, IDN_w~rec_4, IDN_w~rec_5, IDN_w~rec_6, ... 
min values  :          43,          59,          54,          42,          21,           0, ... 
max values  :         787,         665,         641,         543,         640,         742, ...

2. Get data: elevation

We will the function elevation_30s() from geodata. This data comes from the Shuttle Radar Topographic Mission and are aggregated to 1km resolution.

elev <- elevation_30s(country = "IDN")

elev <- crop(elev, extent)

elev
class       : SpatRaster 
dimensions  : 2040, 5424, 1  (nrow, ncol, nlyr)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : 96, 141.2, -11, 6  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source(s)   : memory
name        : IDN_elv_msk 
min value   :        -125 
max value   :        4613

2. Get data: land-cover

Now, let’s get the landcover data. For this we will use the package MODIStsp.

To download MODIS data through the NASA http server, we need to create a profile at https://urs.earthdata.nasa.gov/home to get a user and password.

Today, I’ll present you an example and provide you with the data already processed, but you can do this at home.

MODIStsp_get_prodnames()
  [1] "Surf_Ref_8Days_500m (M*D09A1)"                   
  [2] "Surf_Ref_Daily_005dg (M*D09CMG)"                 
  [3] "Surf_Ref_Daily_500m (M*D09GA)"                   
  [4] "Surf_Ref_Daily_250m (M*D09GQ)"                   
  [5] "Surf_Ref_8Days_250m (M*D09Q1)"                   
  [6] "Snow_Cov_Daily_500m (M*D10A1)"                   
  [7] "Snow_Cov_8-Day_500m (M*D10_A2)"                  
  [8] "Snow_Cov_Day_0.05Deg (M*D10C1)"                  
  [9] "Snow_Cov_8-Day0.05Deg CMG (M*D10C2)"             
 [10] "Snow_Cov_Month_0.05Deg CMG (M*D10CM)"            
 [11] "Surf_Temp_Daily_005dg (M*D11C1)"                 
 [12] "Surf_Temp_Daily_1Km (M*D11A1)"                   
 [13] "Surf_Temp_8Days_1Km (M*D11A2)"                   
 [14] "Surf_Temp_Daily_GridSin (M*D11B1)"               
 [15] "Surf_Temp_8Days_GridSin (M*D11B2)"               
 [16] "Surf_Temp_Monthly_GridSin (M*D11B3)"             
 [17] "Surf_Temp_8Days_005dg (M*D11C2)"                 
 [18] "Surf_Temp_Monthly_005dg (M*D11C3)"               
 [19] "LST_3band_emissivity_Daily_1km (M*D21A1D)"       
 [20] "LST_3band_emissivity_Daily_1km_night (M*D21A1N)" 
 [21] "LST_3band_emissivity_8day_1km (M*D21A2)"         
 [22] "BRDF_Albedo_ModelPar_Daily_500m (MCD43A1)"       
 [23] "BRDF_Albedo_Quality_Daily_500m (MCD43A2)"        
 [24] "Albedo_Daily_500m (MCD43A3)"                     
 [25] "BRDF_Adj_Refl_Daily_500m (MCD43A4)"              
 [26] "BRDF_Albedo_ModelPar_Daily_005dg (MCD43C1)"      
 [27] "BRDF_Albedo_Quality_Daily_005dg (MCD43C2)"       
 [28] "Albedo_Daily_005dg (MCD43C3)"                    
 [29] "BRDF_Adj_Refl_16Day_005dg (MCD43C4)"             
 [30] "AlbPar_1_B1_Daily_30ArcSec (MCD43D01)"           
 [31] "AlbPar_2_B1_Daily_30ArcSec (MCD43D02)"           
 [32] "AlbPar_3_B1_Daily_30ArcSec (MCD43D03)"           
 [33] "AlbPar_1_B2_Daily_30ArcSec (MCD43D04)"           
 [34] "AlbPar_2_B2_Daily_30ArcSec (MCD43D05)"           
 [35] "AlbPar_3_B2_Daily_30ArcSec (MCD43D06)"           
 [36] "AlbPar_1_B3_Daily_30ArcSec (MCD43D07)"           
 [37] "AlbPar_2_B3_Daily_30ArcSec (MCD43D08)"           
 [38] "AlbPar_3_B3_Daily_30ArcSec (MCD43D09)"           
 [39] "AlbPar_1_B4_Daily_30ArcSec (MCD43D10)"           
 [40] "AlbPar_2_B4_Daily_30ArcSec (MCD43D11)"           
 [41] "AlbPar_3_B4_Daily_30ArcSec (MCD43D12)"           
 [42] "AlbPar_1_B4_Daily_30ArcSec (MCD43D13)"           
 [43] "AlbPar_2_B4_Daily_30ArcSec (MCD43D14)"           
 [44] "AlbPar_3_B4_Daily_30ArcSec (MCD43D15)"           
 [45] "AlbPar_1_B5_Daily_30ArcSec (MCD43D16)"           
 [46] "AlbPar_2_B5_Daily_30ArcSec (MCD43D17)"           
 [47] "AlbPar_3_B5_Daily_30ArcSec (MCD43D18)"           
 [48] "AlbPar_1_B6_Daily_30ArcSec (MCD43D19)"           
 [49] "AlbPar_2_B6_Daily_30ArcSec (MCD43D20)"           
 [50] "AlbPar_3_B6_Daily_30ArcSec (MCD43D21)"           
 [51] "AlbPar_1_Vis_Daily_30ArcSec (MCD43D22)"          
 [52] "AlbPar_2_Vis_Daily_30ArcSec (MCD43D23)"          
 [53] "AlbPar_3_Vis_Daily_30ArcSec (MCD43D24)"          
 [54] "AlbPar_1_NIR_Daily_30ArcSec (MCD43D25)"          
 [55] "AlbPar_2_NIR_Daily_30ArcSec (MCD43D26)"          
 [56] "AlbPar_3_NIR_Daily_30ArcSec (MCD43D27)"          
 [57] "AlbPar_1_SWIR_Daily_30ArcSec (MCD43D28)"         
 [58] "AlbPar_2_SWIR_Daily_30ArcSec (MCD43D29)"         
 [59] "AlbPar_3_SWIR_Daily_30ArcSec (MCD43D30)"         
 [60] "BRDF_Albedo_Quality_Daily_30ArcSec (MCD43D31)"   
 [61] "BRDF_Albedo_SolNoon_Daily_30ArcSec (MCD43D32)"   
 [62] "Alb_ValObs_B1_Daily_30ArcSec (MCD43D33)"         
 [63] "Alb_ValObs_B2_Daily_30ArcSec (MCD43D34)"         
 [64] "Alb_ValObs_B3_Daily_30ArcSec (MCD43D35)"         
 [65] "Alb_ValObs_B4_Daily_30ArcSec (MCD43D36)"         
 [66] "Alb_ValObs_B5_Daily_30ArcSec (MCD43D37)"         
 [67] "Alb_ValObs_B6_Daily_30ArcSec (MCD43D38)"         
 [68] "Alb_ValObs_B7_Daily_30ArcSec (MCD43D39)"         
 [69] "BRDF_Albedo_Snow_Daily_30ArcSec (MCD43D40)"      
 [70] "BRDF_Alb_Unc_Daily_30ArcSec (MCD43D41)"          
 [71] "BRDF_Alb_BSA_B1_Daily_30ArcSec (MCD43D42)"       
 [72] "BRDF_Alb_BSA_B2_Daily_30ArcSec (MCD43D43)"       
 [73] "BRDF_Alb_BSA_B3_Daily_30ArcSec (MCD43D44)"       
 [74] "BRDF_Alb_BSA_B4_Daily_30ArcSec (MCD43D45)"       
 [75] "BRDF_Alb_BSA_B5_Daily_30ArcSec (MCD43D46)"       
 [76] "BRDF_Alb_BSA_B6_Daily_30ArcSec (MCD43D47)"       
 [77] "BRDF_Alb_BSA_B7_Daily_30ArcSec (MCD43D48)"       
 [78] "BRDF_Alb_BSA_Vis_Daily_30ArcSec (MCD43D49)"      
 [79] "BRDF_Alb_BSA_NIR_Daily_30ArcSec (MCD43D50)"      
 [80] "BRDF_Alb_BSA_SWIR_Daily_30ArcSec (MCD43D51)"     
 [81] "BRDF_Alb_WSA_B1_Daily_30ArcSec (MCD43D52)"       
 [82] "BRDF_Alb_WSA_B2_Daily_30ArcSec (MCD43D53)"       
 [83] "BRDF_Alb_WSA_B3_Daily_30ArcSec (MCD43D54)"       
 [84] "BRDF_Alb_WSA_B4_Daily_30ArcSec (MCD43D55)"       
 [85] "BRDF_Alb_WSA_B5_Daily_30ArcSec (MCD43D56)"       
 [86] "BRDF_Alb_WSA_B6_Daily_30ArcSec (MCD43D57)"       
 [87] "BRDF_Alb_WSA_B7_Daily_30ArcSec (MCD43D58)"       
 [88] "BRDF_Alb_WSA_Vis_Daily_30ArcSec (MCD43D59)"      
 [89] "BRDF_Alb_WSA_NIR_Daily_30ArcSec (MCD43D60)"      
 [90] "BRDF_Alb_WSA_SWIR_Daily_30ArcSec (MCD43D61)"     
 [91] "BRDF_Albedo_NBAR_Band1_Daily_30ArcSec (MCD43D62)"
 [92] "BRDF_Albedo_NBAR_Band2_Daily_30ArcSec (MCD43D63)"
 [93] "BRDF_Albedo_NBAR_Band3_Daily_30ArcSec (MCD43D64)"
 [94] "BRDF_Albedo_NBAR_Band4_Daily_30ArcSec (MCD43D65)"
 [95] "BRDF_Albedo_NBAR_Band5_Daily_30ArcSec (MCD43D66)"
 [96] "BRDF_Albedo_NBAR_Band6_Daily_30ArcSec (MCD43D67)"
 [97] "BRDF_Albedo_NBAR_Band7_Daily_30ArcSec (MCD43D68)"
 [98] "Vegetation_Indexes_16Days_500m (M*D13A1)"        
 [99] "Vegetation_Indexes_16Days_1Km (M*D13A2)"         
[100] "Vegetation_Indexes_Monthly_1Km (M*D13A3)"        
[101] "Vegetation_Indexes_16Days_005dg (M*D13C1)"       
[102] "Vegetation_Indexes_Monthly_005dg (M*D13C2)"      
[103] "Vegetation Indexes_16Days_250m (M*D13Q1)"        
[104] "LAI_8Days_500m (MCD15A2H)"                       
[105] "LAI_4Days_500m (MCD15A3H)"                       
[106] "LAI_8Days_500m (M*D15A2H)"                       
[107] "Net_ET_8Day_500m (M*D16A2)"                      
[108] "Net_ETgf_8Day_500m (M*D16A2GF)"                  
[109] "Net_ETgf_Yearly_500m (M*D16A3GF)"                
[110] "Gross_PP_8Days_500m (M*D17A2H)"                  
[111] "Gross_PP_GapFil_8Days_500m (M*D17A2HGF)"         
[112] "Net_PP_GapFil_Yearly_500m (M*D17A3HGF)"          
[113] "Veg_Cont_Fields_Yearly_250m (MOD44B)"            
[114] "Land_Wat_Mask_Yearly_250m (MOD44W)"              
[115] "Burned_Monthly_500m (MCD64A1)"                   
[116] "ThermalAn_Fire_Daily_1Km (M*D14A1)"              
[117] "ThermalAn_Fire_8Days_1Km (M*D14A2)"              
[118] "LandCover_Type_Yearly_005dg (MCD12C1)"           
[119] "LandCover_Type_Yearly_500m (MCD12Q1)"            
[120] "LandCover_Dynamics_Yearly_500m (MCD12Q2)"        
[121] "Dwnwrd_Srw_Rad_3h_005dg (MCD18A1)"               
[122] "Dwnwrd_PAR_3h_005dg (MCD18A2)"                   
[123] "MAIA_Land_Surf_BRF (MCD19A1)"                    
[124] "MAIA_Land_AOT_daily (MCD19A2)"

We will use the Land Cover Type Yearly L3 Global 500m, and download data for Indonesia from the year 2020.

2. Get data: land-cover

With MODIStsp_get_prodlayers() you can see all the layers of the product:

MODIStsp_get_prodlayers("MCD12Q1", prodver = "061")$bandnames
 [1] "LC1"                 "LC2"                 "LC3"                
 [4] "LC4"                 "LC5"                 "LC_Prop1"           
 [7] "LC_Prop2"            "LC_Prop3"            "LC_Prop1_Assessment"
[10] "LC_Prop2_Assessment" "LC_Prop3_Assessment" "LC_QC"              
[13] "LC_LW"


There are five different land cover classification schemes, we will be using the primary land cover scheme (LC1) which identifies 17 classes defined by the IGBP (International Geosphere-Biosphere Programme), including 11 natural vegetation classes, 3 human-altered classes, and 3 non-vegetated classes.

2. Get data: land-cover

We will use the function MODIStsp() to download data:

MODIStsp(
  ...,
  gui = TRUE,
  out_folder = NULL,
  out_folder_mod = NULL,
  opts_file = NULL,
  selprod = NULL,
  prod_version = NULL,
  bandsel = NULL,
  quality_bandsel = NULL,
  indexes_bandsel = NULL,
  sensor = NULL,
  download_server = NULL,
  downloader = NULL,
  user = NULL,
  password = NULL,
  download_range = NULL,
  start_date = NULL,
  end_date = NULL,
  spatmeth = NULL,
  start_x = NULL,
  end_x = NULL,
  start_y = NULL,
  end_y = NULL,
  bbox = NULL,
  spafile = NULL,
  out_projsel = NULL,
  output_proj = NULL,
  out_res_sel = NULL,
  out_res = NULL,
  resampling = NULL,
  reprocess = NULL,
  delete_hdf = NULL,
  nodata_change = NULL,
  scale_val = NULL,
  ts_format = NULL,
  out_format = NULL,
  compress = NULL,
  test = NULL,
  n_retries = 5,
  verbose = TRUE,
  parallel = TRUE
)

2. Get data: land-cover

We used the function MODIStsp() to download data for you. MODIS require to create a user and password.

#| echo: true
#| eval: false
MODIStsp(
  gui = FALSE,
  out_folder = "data/",
  out_folder_mod = "data/",
  selprod = "LandCover_Type_Yearly_500m (MCD12Q1)",
  prod_version = "061",
  bandsel = "LC1",
  sensor = "Terra",
  user = your_user, # your username for NASA http server
  password = your_pass, # your password for NASA http server
  start_date = "2020.01.01",
  end_date = "2020.12.31",
  verbose = TRUE,
  bbox = extent, # bbox covering Indonesia
  spatmeth = "bbox",
  out_format = "GTiff",
  compress = "LZW",
  out_projsel = "User Defined",
  output_proj = 4326,
  delete_hdf = TRUE,
  parallel = TRUE
)

2. Get data: land-cover

Let’s read the processed layer.

land <- rast("data/LandCover_Type_Yearly_500m_v61/LC1/MCD12Q1_LC1_2020_001.tif")

land <- crop(land , extent)

land
class       : SpatRaster 
dimensions  : 4244, 11483, 1  (nrow, ncol, nlyr)
resolution  : 0.004005922, 0.004005655  (x, y)
extent      : 96, 142, -11, 6  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : MCD12Q1_LC1_2020_001.tif 
name        : MCD12Q1_LC1_2020_001

2. Get data: anthropogenic

Again using the package geodata.

footp <- footprint(year = 2009)

footp <- crop(footp, extent)

footp
class       : SpatRaster 
dimensions  : 2040, 5520, 1  (nrow, ncol, nlyr)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : 96, 142, -11, 6  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source(s)   : memory
name        : wildareas-v3-2009-human-footprint 
min value   :                                 0 
max value   :                                49

Let’s compare our rasters

Rasters should have the same spatial extent and resolution…

If they don’t match… what should we do?

compareGeom(tavg, prec, elev, footp, land, stopOnError = F, messages = T, res = T)
[1] FALSE

Let’s compare our rasters

Rasters should have the same spatial extent and resolution…

If they don’t match… what should we do?

st_bbox(tavg)
 xmin  ymin  xmax  ymax 
 96.0 -11.0 141.5   6.0 
st_bbox(prec) 
 xmin  ymin  xmax  ymax 
 96.0 -11.0 141.5   6.0 
st_bbox(elev) 
 xmin  ymin  xmax  ymax 
 96.0 -11.0 141.2   6.0 
st_bbox(footp)
xmin ymin xmax ymax 
  96  -11  142    6 
st_bbox(land)
xmin ymin xmax ymax 
  96  -11  142    6 
#Compare with the country extent
st_bbox(country)
      xmin       ymin       xmax       ymax 
 95.293026 -10.359987 141.033852   5.479821

2.1. Process raster

2.1. Process raster

We will use the package terra to process the rasters.
First we need to create a new raster to use as a template.

r <- rast(res = 10000, # Target cell size
          ext(vect(st_transform(country, project_crs))), # Extent of the raster
          crs = project_crs) %>% # Useful planar projection
  project(., "epsg:4326") # Projection of our rasters
r
class       : SpatRaster 
dimensions  : 176, 508, 1  (nrow, ncol, nlyr)
resolution  : 0.08993314, 0.08993314  (x, y)
extent      : 95.29303, 140.9791, -10.38937, 5.438861  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)

2.1. Process raster

We will resample() all of the rasters to our target resolution.

tavg <- resample(tavg, r, method = "bilinear")

|---------|---------|---------|---------|
=========================================
                                          
prec <- resample(prec, r, method = "bilinear")

|---------|---------|---------|---------|
=========================================
                                          
elev <- resample(elev, r, method = "bilinear")
footp <- resample(footp, r, method = "near")
land <- resample(land, r, method = "near")

And we can compare the rasters after resampling.

compareGeom(tavg, prec, elev, footp, land, stopOnError = F, messages = T, res = T)
[1] TRUE

2.1. Process raster

And finally we will use mask() to intersect the layer with Indonesia.

tavg <- mask(tavg, vect(country))
prec <- mask(prec, vect(country))
elev <- mask(elev, vect(country))
footp <- mask(footp, vect(country))
land <- mask(land, vect(country))

2.2. Visualise data

2.2. Visualise data

We will use the package tmap.

tmap_mode(mode = "view")

2.2. Visualise data

Annual mean temperature (WorldClim BIO1)

tm_shape(mean(tavg)) +
  tm_raster(palette = "OrRd", midpoint = NA, style = "cont", n = 10)

2.2. Visualise data

Precipitation seasonality (WorldClim BIO15)

tm_shape(mean(prec)) +
  tm_raster(palette = "RdBu", midpoint = NA, style = "cont", n = 10) +
  tm_view(bbox = extent)

2.2. Visualise data

Elevation (WorldClim)

tm_shape(elev) +
  tm_raster(palette = "BrBG", midpoint = NA, style = "cont", n = 10) +
  tm_view(bbox = extent)

2.2. Visualise data

Land-cover (MODIS Terra MCD12Q1). We need to make sure values are factors.

names(land)
[1] "MCD12Q1_LC1_2020_001"
land <- as.factor(land$MCD12Q1_LC1_2020_001)

Then we can add values to the levels.

levels(land[[1]]) <- data.frame(
  ID = c(1, 2, 4, 5, 6, 8, 9, 10, 11, 12, 13, 14, 16, 17),
  label = c(
    "Evergreen Needleleaf Forest",
    "Evergreen Broadleaf Forests",
    "Deciduous Broadleaf Forests",
    "Mixed Forests",
    "Closed Shrublands",
    "Woody Savannas",
    "Savannas",
    "Grassslands",
    "Permanent Wetlands",
    "Croplands",
    "Urban and Built-up Lands",
    "Cropland/Natural Vegetation Mosaics",
    "Barren",
    "Water Bodies"
  )
)

2.2. Visualise data

Land-cover (MODIS Terra MCD12Q1)

tm_shape(land) +
  tm_raster(palette = "Set1", style = "cat") +
  tm_view(bbox = extent)

2.2. Visualise data

Human footprint index (SEDAC)

tm_shape(footp) +
  tm_raster(palette = "viridis", style = "cont") +
  tm_view(bbox = extent)

3. Correlations between predictors

3. Correlations between predictors

Let’s get values per grid-cell.

names(tavg)
 [1] "IDN_wc2.1_30s_tavg_1"  "IDN_wc2.1_30s_tavg_2"  "IDN_wc2.1_30s_tavg_3" 
 [4] "IDN_wc2.1_30s_tavg_4"  "IDN_wc2.1_30s_tavg_5"  "IDN_wc2.1_30s_tavg_6" 
 [7] "IDN_wc2.1_30s_tavg_7"  "IDN_wc2.1_30s_tavg_8"  "IDN_wc2.1_30s_tavg_9" 
[10] "IDN_wc2.1_30s_tavg_10" "IDN_wc2.1_30s_tavg_11" "IDN_wc2.1_30s_tavg_12"

3. Correlations between predictors

Let’s see all the values together.

predictors <- tibble(
  "tavg" = values(tavg$IDN_wc2.1_30s_tavg_1) %>% as.vector(),
  "prec" = values(prec$IDN_wc2.1_30s_prec_1) %>% as.vector(),
  "elev" = values(elev) %>% as.vector(),
  "land" = values(land) %>% as.vector(),
  "footp" = values(footp) %>% as.vector()
) %>% 
  mutate(
    across(everything(), ~replace_na(.x, NA))
  )

3. Correlations between predictors

Let’s see all the values together.

summary(predictors)
      tavg            prec             elev              land      
 Min.   : 7.02   Min.   : 86.85   Min.   :   1.14   Min.   : 1.0   
 1st Qu.:24.31   1st Qu.:236.56   1st Qu.:  34.65   1st Qu.: 2.0   
 Median :25.90   Median :283.48   Median : 131.18   Median : 2.0   
 Mean   :25.00   Mean   :278.73   Mean   : 355.91   Mean   : 6.1   
 3rd Qu.:26.71   3rd Qu.:316.51   3rd Qu.: 479.28   3rd Qu.: 9.0   
 Max.   :28.48   Max.   :680.45   Max.   :3931.85   Max.   :17.0   
 NA's   :71703   NA's   :71703    NA's   :71849     NA's   :69753  
     footp      
 Min.   : 0.00  
 1st Qu.: 1.00  
 Median : 5.00  
 Mean   : 6.76  
 3rd Qu.:11.00  
 Max.   :47.00  
 NA's   :69753

3. Correlations between predictors

Let’s see how these values correlate.

What can this plot tell you?

EXTRA Can we explore any pattern with the hotspots of turtle species richness from last class?

Turtles of Asia

Let’s read the data, plot it and see the patterns

grid_testudines_Asia <- readRDS("data/grid_testudines_Asia.rds") %>%
  st_transform(., 4326) %>%
  st_crop(country)

Turtles of Asia

Let’s read the data, plot it and see the patterns

tm_shape(grid_testudines_Asia) +
  tm_fill("SR", palette = "YlOrBr") +
  tm_view(bbox = extent)

Turtles of Asia

Let’s read the data, plot it and see the patterns

Turtles of Asia

Let’s read the data, plot it and see the patterns

Any doubts?